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Abstract 

This paper presents a data-driven receding horizon fault estimation method for additive actuator and sensor faults in unknown 
linear time-invariant systems, with enhanced robustness to stochastic identification errors. State-of-the-art methods construct 
fault estimators with identified state-space models or Markov parameters, but they do not compensate for identification 
errors. Motivated by this limitation, we first propose a receding horizon fault estimator parameterized by predictor Markov 
parameters. This estimator provides (asymptotically) unbiased fault estimates as long as the subsystem from faults to outputs 
has no unstable transmission zeros. When the identified Markov parameters are used to construct the above fault estimator, 
zero-mean stochastic identification errors appear as model uncertainty multiplied with unknown fault signals and online system 
inputs/outputs (I/O). Based on this fault estimation error analysis, we formulate a mixed-norm problem for the offline robust 
design that regards online I/O data as unknown. An alternative online mixed-norm problem is also proposed that can further 
reduce estimation errors when the online I/O data have large amplitudes, at the cost of increased computational burden. 
Based on a geometrical interpretation of the two proposed mixed-norm problems, systematic methods to tune the user-defined 
parameters therein are given to achieve desired performance trade-offs. Simulation examples illustrate the benefits of our 
proposed methods compared to recent literature. 
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1 Introduction 

Model-based fault diagnosis techniques for linear dy¬ 
namic systems have been well established during the past 
two decades [2,5,7,16]. Recently, the model-based re¬ 
ceding horizon approach has received attention because 
it provides a flexible framework to enhance robustness 
of passive fault diagnosis [34,36] and to enable optimal 
input design in active fault diagnosis [28,29,31]. How¬ 
ever, an explicit and accurate system model is often un¬ 
known in practice. In such situations, a conventional ap¬ 
proach first identifies the system model from system I/O 
data, and then designs the model-based fault diagnosis 
system under various performance criteria [22,25,32]. 
Without explicitly identifying a system model, recent 
research efforts investigate data-driven approaches to 
construct a fault diagnosis system utilizing the link be- 
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tween system identification and the model-based fault 
diagnosis methods [8,9,30]. These recent data-driven 
approaches simplify the design procedure by skipping 
the realization of an explicit system model, while at the 
same time allow developing systematic methods to ad¬ 
dress the same fault diagnosis performance criteria as 
the existing model-based approaches. 

Most recent data-driven fault diagnosis approaches 
for unknown linear dynamic systems can be classi¬ 
fied into two categories. The first category, e.g., [27] 
and [9, 10], identifies a projection matrix known as 
parity space/vectors for residual generation, by ex¬ 
ploiting the subspace identification method based on 
principal component analysis (SIM-PCA) [17]. How¬ 
ever, as pointed out in [13], a model reduction step is 
needed to determine the projection matrix, hence leads 
to the nonlinear dependence of the generated residuals 
on the identification errors. Therefore it is difficult to 
guarantee the robustness of such data-driven methods 
to the identification errors. 

The second category of data-driven fault diagnosis 
methods, e.g., [11], utilizes the Markov parameters (or 
impulse response parameters) which can be obtained in 


Preprint submitted to Automatica 


2 March 2015 



the first step of the predictor based subspace identifi¬ 
cation (PBSID) technique [6,33]. It constructs residual 
generators parameterized by the predictor Markov pa¬ 
rameters. The main advantage of this method is that 
the residual signal linearly depends on the identification 
errors of the predictor Markov parameters. Hence a ro¬ 
bust scheme has been developed in [13,14] to cope with 
stochastic identification errors. This benefit of robust¬ 
ness compared to the SIM-PCA based method in [10] is 
achieved at the cost of increased computational burden 
in incorporating past I/O data. 

Most of the data-driven fault diagnosis literature men¬ 
tioned above discuss only fault detection and isolation. 
It is much more involved to estimate/identify the fault 
signal in the data-driven setting. The work in [1, 26] 
proposed to reconstruct faults by minimizing the recon¬ 
structed squared prediction error obtained from PCA. 
However, this approach did not fully investigate the sta¬ 
tistical properties of the calculated fault estimates. By 
investigating the link between system-inversion based 
fault reconstruction and the predictor Markov param¬ 
eters, the method in [12] constructed fault estimators 
parameterized by the predictor Markov parameters. Its 
fault estimates are asymptotically unbiased as the esti¬ 
mation horizon length tends to infinity, under the con¬ 
dition that the underlying inverted system is stable. 

One drawback of the data-driven fault estimator pro¬ 
posed in [12] is that it cannot be directly applied to 
sensor faults in an unstable open-loop plant because its 
underlying inverted system is unstable. Another limita¬ 
tion of this method is that it does not compensate for the 
identification errors. The robustness of fault estimation 
to the identification errors is critical in two situations: 
1) there exist large identification errors due to small 
number of identification data samples or low signal-to- 
noise ratio in identification data; 2) multiplication of the 
erroneous identified matrices with online I/O data of 
large amplitude cannot be simply ignored. 

Motivated by the above two drawbacks of the proposed 
method in [12], this paper develops data-driven robust 
fault estimation methods for additive actuator/sensor 
faults, utilizing the identified Markov parameters. In 
order to pave the way for data-driven design, we first 
construct a receding horizon (RH) fault estimator pa¬ 
rameterized by the predictor Markov parameters, as¬ 
suming that the predictor Markov parameters are accu¬ 
rately available. It gives (asymptotically) unbiased fault 
estimates under the condition that the subsystem from 
faults to outputs has no unstable transmission zeros. The 
above condition for unbiasedness generalizes the require¬ 
ment of stable inversion in [12]. An immediate benefit is 
that our fault estimator can be applied to sensor faults 
in unstable open-loop plants as long as the above condi¬ 
tion for unbiasedness is satisfied, whereas the proposed 
method in [12] cannot. 


Our data-driven design parameterizes the above RH 
fault estimator with predictor Markov parameters iden¬ 
tified from closed-loop data. The obtained data-driven 
fault estimation error is linear with regards to the sto¬ 
chastic identification errors of Markov parameters, al¬ 
though the identification errors appear as multiplicative 
uncertainty that couples with unknown fault signals 
as well as online I/O data. In order to enhance ro¬ 
bustness to stochastic identification errors, we propose 
two mixed-norm fault estimators. The first one can be 
designed offline by regarding the online I/O data as un¬ 
known. By exploiting online I/O data in its formulated 
mixed-norm problem, the second robust fault estimator 
further reduces estimation errors when the online I/O 
data have large amplitudes, at the cost of increased 
online computational burden. Based on a geometric in¬ 
terpretation of the formulated mixed-norm problems, a 
systematic tuning method for the user-defined parame¬ 
ters therein is provided to achieve the desired trade-offs 
between estimation bias and variance. Our proposed 
methods can handle sensor and actuator faults either 
separately or simultaneously. Only the separate scenario 
is illustrated in detail in this paper. Exact formulas for 
the simultaneous scenario can be derived in a straight¬ 
forward manner but are omitted for the sake of brevity. 

The rest of this paper starts with the problem formula¬ 
tion and some preliminaries on closed-loop identification 
of predictor Markov parameters in Section 2. Section 3 
constructs the predictor-based RH fault estimator, and 
analyzes its condition for unbiasedness. A data-driven 
nominal fault estimator is given in Section 4. Section 
5 and 6 propose two mixed-norm fault estimators with 
enhanced robustness to identification errors. Simulation 
studies are finally given in Section 7. 


2 Preliminaries and problem formulation 

2.1 Notations 

For a matrix X, its range and null space is denoted by 
TZ(X) and Af (X), respectively. X~ represents the left 
inverse satisfying X~X = /, while X^ represents the 
generalized inverse satisfying 

XX (1) X = X. (1) 

X W represents the i th column of X. The trace of X is 
denoted by tr(X). Let ||X|| F represent the Frobenius 
norm of the matrix X. The minimal eigenvalue of a sym¬ 
metric matrix X is represented by A m ; n (X). Let vec (X) 
represent the column vector concatenating the columns 
of a matrix X. The symbol “(g)” stands for Kronecker 
product. Let diag (Xi, X 2 , • • • , X n ) denote a block di¬ 
agonal matrix with Xi, X 2 , • • • , X n as its diagonal ma¬ 
trices. 
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2.2 Problem formulation 

We consider linear discrete-time systems governed by 
the following state space model: 

£(fc + 1) = A£(k) + Bu(k) + Ef(k) + Fw(k) 

y(k) = c{(k) + Du(k) + Gf(k) + v(k). } 

Here £(/c) G R n , y(k) G and u(k) G R" u rep¬ 

resent the state, the output measurement, and the 
known control input at time instant fc, respectively. 
The process and measurement noises w(k) G M"” and 
v[k) G K"’' are white zero-mean Gaussian, with covari¬ 
ance matrices E (w(k)w T (k)) = Q, E (v(k)v T (k)) = R, 
E (w(k)v T (fc)) = 0. f(k) G is the unknown fault 
signal to be estimated. A, B, C, D, E, F, G are constant 
real matrices, with bounded norms and appropriate 
dimensions. 

The following assumption is standard in Kalman fil¬ 
tering [18] and subspace identification [6,19]: 

Assumption 1 The pair (C, A) is assumed detectable; 
and there are no uncontrollable modes of ^A, FQi'j on 

the unit circle, where Q s • (Q^ = Q is the covariance 

matrix of w(k). 

Based on Assumption 1, the system (2) admits the one- 
step-ahead predictor form given by [18] 

x(k + 1) = $x(fc) + Bu(k) + Ef(k) + Ky(k) , . 
y(k) = Cx(k ) + Du(k) + Gf(k) + e(k), 

where K is the steady-state Kalman gain, $ = A — KC, 
B = B — KD , and E = E — KG, |e(fc)} is the zero-mean 
innovation process with the covariance matrix £ e . 

We consider additive sensor or actuator faults in this 
paper, i.e., 

• fault of the j th sensor: 


E = 0 naXl ,G = lW, E = -KW; (4) 

• fault of the I th actuator: 

E = B [l] , G = D [l \ E = B [l] ; (5) 

• simultaneous faults of the j th sensor and I th actuator: 


E = 


On, XI BM 


,G = 


j\j] r>v\ 


,E = 


~kW bN 

( 6 ) 


with X representing the j th column of a matrix X. 
Denote the predictor Markov parameters by 


m 


o 

II 

f 0 

i = 0 

,H¥ = < 

) 

\ 


C^-'B i > 0 

{ C<S> 1 ~ 

i > 0 



G i = 0 
C&^Ei >0 


( 7 ) 


Assumption 2 The relative degree of the fault sub¬ 
system E, C, G^J is t, i.e., r is the smallest nonneg- 

f f f 

ative integer i such that Hq = H[ = • • • = H(_ 1 = 0 
and H- ^ 0 [20]; moreover, rank (Hf) = nf [12]. 


Note that r = 0 for sensor faults and r > 0 for actuator 
faults. 


The essential goals of this paper are to design a fault 
estimator from identification data without knowing the 
system matrices in (2), and moreover to robustify the 
fault estimator against identification errors. 

Concerning the identification data, it should be noted 
that in practice data from faulty conditions may be sel- 
domly available, or if recorded then without a reliable 
fault description [9]. Hence we make the assumption as 
below: 

Assumption 3 Only I/O data collected from the fault- 
free condition are used in our data-driven design. 

In contrast to [24] which assumes the fault signals f(k) 
evolve according to a random walk model, no assumption 
is made in this paper about how the fault signals f(k) 
vary with time. 

2.3 Closed-loop identification of predictor Markov pa¬ 
rameters 


Considering Assumption 3, we set f(k) = 0 in (2) for 
the identification data collected from the fault-free con¬ 
dition. Then with f(k) = 0, the predictor form (3) over 
the time window [t, ■ ■ ■ , t + N — 1] can be written into 
the following data equation [6,33]: 


Y id = C'$ p X id + SZ id + E id , (8) 


where 


TTU Try 
**p 


H f H\ H/ 


(9) 
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denotes the sequence of Markov parameters {H/} and 
{Hf} (defined in (7)) to be identified. The detailed def¬ 
initions of the data matrices X; d , Yj d and Zj d can be 



found in [33], and Ej d is the sequence of the innovation 
signal in the identification data. 

The least-squares (LS) estimate of the Markov parame¬ 
ters S is 


with ko = k — L + 1. For the predictor form (3), let Ol 
denote its extended observability matrix with L block 
elements, and be the lower triangular block-Toeplitz 

matrix with L block columns and rows, with * repre¬ 
senting u, y , or /: 


2 = arg_min ||Y id - SZ id ||^ = Y id Z id 
= Z + C$ p X id Z7 d + E id Zy, 


( 10 ) 


with Ztj = (Z id Z^i) . As standard assumptions 

for consistent identification from closed-loop data, we 
assume that 1) the data matrix Z; d has full row rank, and 
2) either the controller has at least one-step delay or the 
plant model has no direct feedthrough (D = 0) [6,33]. 


With sufficiently large p, the estimation bias C ,< I» J 'Xi d Z^ 
can be neglected. Then the stochastic identification er¬ 
rors are 

AS = S —H«E id Zy. (11) 

Hence according to (11), the identification errors in 
Markov parameters can also be written as 


A H? = H? - H? = E id M“, 
A HV = Hf- H\ - E id Mf, 


( 12 ) 


c 


’ HI 

o . 

.. 0 

CA> 

H 

ti* 

II 

m 

... tnj 

O* 

0 



1 

HI- 2 • • 

• H o. 


(16) 

Given the I/O data over the sliding window [ko, k ], the 
stacked residual signal r k L in [ko, k\ can be computed by 

rfc,L = yk,L — T^y fc)i — T£ufc,L, (17) 

according to the predictor form (3). We can further write 
down the transitions from unknown initial state, faults 
and noises to the stacked residual signal r kL as 


r k,L = O L x(k 0 ) + T{f fc , L + e fc;L . (18) 


where Hf and Hf represent the estimated Markov pa¬ 
rameters in H given by (10), M“ and M? are the corre¬ 
sponding blocks of Z ;d , i.e., 


Z 


id 


M™ M« ■■■ M“ Mf Mtf 


, Ml = 0. (13) 


The innovation covariance can be estimated by [16,19] 

E e = cov (Y id - SZ id ) . (14) 

For the sake of brevity, we shall not distinguish between 
the estimated innovation covariance E e and its true value 
E e in the rest of this paper. 


With Assumption 2, (18) can be simplified as 


r k,L 


Ol t{ >t 

1 

5 

i_ 


I I/c— t,L—t \ 


*L,t v -v~ 


+ e fc,Lj 


(19) 


where r is the relative degree of the fault subsystem 
(A,E,C,G), T j l T represents the first L — r block- 

columns of defined similar to (16), f k - T ,L-T is 
defined in the same way as in (15). 

With (19), we can formulate the receding horizon fault 
estimation (RHFE) problem 


3 Predictor-based receding horizon fault esti¬ 
mation 

In this section, we will construct an RH fault estimator 
based on the predictor form of the system (2). Here we 
consider the predictor form instead of the original system 
model (2) in order to pave the way for data-driven design. 


min ||r feiL - T L , r f/ : _ Tii _ r ||^_ 1 (20) 

k — t , L — t e,L 

in the LS sense, with 

E e ,i = II <S> E e (21) 


Consider a sliding window with a length of L sampling 
instants. Define stacked data vectors in this window as 
u k,L, y k,L, f k,L, and e k ,L, respectively for the signals u, 
V , f, and e; e.g., 


denoting the covariance matrix of e k ,L- It has non-unique 
solutions because ^ l,t may not have full column rank. 
One solution to the problem (20) is 


u k,L = 


(k 0 ) 


C k ) 


(15) 


L k—r,L—r 


= * 




L ,t e,L 


(i) 


* 


L,rS"ir kt L. (22) 
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We will show in the following theorem, however, that the 
last rif entries of fj?_ T L _ T , he., 

f(k- T )=l nf fLr,L-r (23) 

with I nf = [0 I nf ] € K n / x ( n + n /( L —r)), represent an 

(asymptotically) unbiased estimate of / (k — r) under 
certain conditions. The estimation delay r in (23) is 
caused by the relative degree in Assumption 2. 

Theorem 1 Let r and v denote the relative degree 
and the observability index of the fault subsystem 
(&,E,C,G), respectively. 


with the unknown input. Although it seems that the RHFE 
problem (20) jointly estimates initial state and faults, we 
are actually only interested in the fault estimate without 
unbiased reconstruction of the unknown initial state. This 
is an intuitive reason why Theorem 1 can cope with the 
unknown initial state in the case that ($, C) is detectable. 

Remark 2 Theorem 1 above generalizes Theorems 1 
and 2 in [35] in two aspects: 1) Theorems 1 and 2 in [35] 
are limited to the case t = 0, while Theorem 1 here 
applies to general relative degrees; 2) Theorems 1 and 2 
in [35] focus on the fault estimator constructed with the 
original system (2), while in this work we construct in 
Theorem 1 the fault estimator with the predictor (3). 


(i) The r-delay fault estimate f(k — r) defined in 
(23) is unbiased for all L > v + r if and only if 
(<h, E , O t+ i, H [) has no transmission zeros, with 


It should be noted that an RHFE problem similar to 
(20) can also be formulated using the original system 
(2), see [35]. Its equivalence to our RHFE problem (20) 
is shown in the following theorem. 


H 1 = 


wr {h{? 


( H ]) 1 


lT 


(24) 


(ii) The t- delay fault estimate f(k—r) is asymptotically 
unbiased for L —> oo if and only if all transmission 
zeros of ( ( f>, E, O t+ i, H([) are stable. 


Theorem 2 If both the original system model (2) and its 
predictor form (3) are accurately available, the r-delay 
fault estimate f(k — r), computed by (22) and (23) based 
on the predictor form (3), is equivalent to the fault es¬ 
timate proposed as Equation (15) in [35] based on the 
original system model (2). 


The proof is given in Appendix B. 

Instead of including the unknown initial state as in the 
RHFE problem (20), the essential idea of [12] is to find 
a lower triangular block-Toeplitz matrix Tf such that 
T 9 l ■ T[ = I and the estimation error caused by the 
unknown initial state exponentially decays with L. The 
condition for unbiasedness in [12] requires that the in¬ 
verse system related to T® is stable. However this has 
several drawbacks: it does not clarify how the unbiased¬ 
ness condition is related to the system property of the 
underlying plant; and moreover, for the case of sensor 
faults in an open-loop unstable plant, the method in [12] 
cannot find a stable left inverse matrix T 9 L for T^ . 

On the contrary, Theorem 1 clearly states that the con¬ 
dition for unbiasedness is related to the invariant zeros 
of the fault subsystem in the underlying plant. An imme¬ 
diate benefit is that our proposed RH fault estimator can 
ensure (asymptotically) unbiased estimates for sensor 
faults in an open-loop unstable plant, as long as the fault 
subsystem has no unstable transmission zeros. 

Remark 1 The unbiasedness condition of the r-delay 
fault estimate stated in Theorem 1 has close links with 
the r-delay left, inversion in [15,23] and the r-delay input 
and initial-state reconstruction in [20]. However, the r- 
delay left inversion in [15, 23] requires the initial state 
to be known a priori, while the r-delay input and initial- 
state reconstruction in [20] requires observability of the 
pair ($, C) to simultaneously reconstruct the initial state 


The proof of Theorem 2 is given in Appendix C. The 
above equivalence implies that the predictor gain K does 
not affect the statistics of the fault estimation error, and 
the condition of unbiasedness in Theorem 1 holds for the 
RH fault estimation using the original form. 


4 Data-driven nominal receding horizon fault 
estimator 


In this section, we will parameterize the RH fault esti¬ 
mator introduced in Section 3 with the predictor Markov 
parameters, and then provide the nominal data-driven 
design method without considering identification errors. 

In order to construct the LS fault estimator (22), we first 
need to construct the block-Toeplitz matrices Tf, T V L , 

f 

and Tf from the predictor Markov parameters according 
to (16). Then, we need the extended observability matrix 
Ol • One possible approach is to identify Ol from the 
block-Hankel matrix 


H 


O 

L,m 


H( 

Hf 


m 

Hf 


HU 


H. 


m+1 


(25) 


TTU TTU 

L h l h l +i 


ffE+m-lJ 


through a model reduction step [33]. But this model re¬ 
duction step would make the fault estimation error de¬ 
pend nonlinearly on the identification errors. In order to 
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avoid this difficulty, we substitute C>Lx(k o) = H ° Lm Cm 
into (19) by exploiting the following property: 

n{o L ) = n(w L J ( 26 ) 


for m > n. Then (19) can be rewritten as 


?k,L = 


H 


L,m i L, r 


Cm 

ffc— t,L— t 


+ e fc,L> (27) 


r L , T 


where Tf consists of the first L —r block-columns of 
defined in (16). By doing so, the fault estimation error 
becomes linear with regards to the identification errors, 
as shown later in (43). Based on (27), an LS problem 
similar to (20) can be formulated, and one solution is 


K 

L k—r,L—r 


(i) 


(T£ T E-iT L>T ) TI jT E :lr k , L . (28) 
Similarly to (23), we obtain the fault estimate 

f{k-r)=l nf fi_ T r _ T , (29) 


with I nf = 


0 In 


l f k — T,L — T'> 

G R n /x(n u -m+ra „(£-r))_ 


Algorithm 1 Data-driven nominal RH fault estimation 

1) Collect identification data from the fault-free con¬ 
dition, and form the data matrices Yi d and Z; d 
with sufficiently large p [33]. 

2) Compute the sequence of Markov parameters 5 
and the innovation covariance £ e via (10) and 
(14); extract the identified Markov parameters 
77“ and Hf from 5 according to (9); and extract 
Hl according to (4)-(7): 

• for j th sensor faults: 

H{ = for i > 0, and H f 0 = i w ; 

(32) 

• or for 7 th actuator faults: 

H{ = (H?) [l] (i > 0). (33) 

3) Select sufficiently large L. Construct the esti¬ 
mates of S e , L in (21), T y L , T“, in (16), U° L rn in 
(25), and T l , t in (27) as S e ,L, T" , T U L , T{, 

and Yl T by using £ e and the identified Markov 
parameters {Hf , Hf , if/}. Form t with the 

first L — t block-columns of T;£. 

4) Compute the nominal fault estimator according 
to (30) and (31). 


Theorem 3 The sufficient and necessary condition for 
unbiased estimation in Theorem 1 applies to the fault 
estimate defined in (28)-(29). 

The proof is given in Appendix D. 

Combining (17), (28), and (29) yields the RH fault esti¬ 
mator as below: 


5 Data-driven robust receding horizon fault es¬ 
timation 

The data-driven nominal design in Algorithm 1 might 
give biased fault estimates due to errors in the identified 
Markov parameters. To address this problem, this sec¬ 
tion proposes an offline robust design which regards the 
online I/O data as unknown in the design stage. 


f(k -t) = G n v k)L = Q n 


I - T y L i -T“ 


<?n = Tn f (Tl iT E"iT L , T ) (1) Tl iT E"i, (31) 



y k,L 


. W.L . 


(30) 


5.1 Data-driven robust design 

Since the Markov parameters related to faults are ex¬ 
tracted from Hf or Hf via (33) or (32), the identifica- 

t 

tion errors of Hf can be expressed as 


where Gn represents the nominal RH fault estimator 
based on the residual signal r k,L- 

Without considering the identification errors, the data- 
driven design of nominal RH fault estimator can now be 
summarized in Algorithm 1. For the sake of brevity, we 
do not list the estimated fault Markov parameters H / 
and their estimation errors for simultaneous sensor and 
actuator faults, because they can be straightforwardly 
derived similarly to (32) and (33). Thus all our proposed 
algorithms in this paper can be directly extended to deal 
with simultaneous sensor and actuator faults. 


A H{ = E id M/, (34) 

where 

(.Mf )^ for faults of the 7 th actuator 

v 1 > J (35) 

— (Mf)^ for faults of the j th sensor 

with Mf and Mf defined in (12)-(13). 

With (12) and (34), the estimated matrices Tf, Tf, 
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TL, H? m and f L , T in Algorithm 1 can be written as 

H lm = H ° Lt m + EidMT* = T“ - E id M*, (36) 
Tl = T l + E id M l, T [ T = T f LT + E id M{ jr , (37) 
Tl, t = T iir + E id M T , (38) 


where M 1 / m is the block-Hankel matrix constructed 
with M“, Mg, • • ■ , Mg +m _ l similarly to H ° Lm hr (25), 
M x is the block-Toeplitz matrix constructed with 
Mg, M*, • • • , M£_ 1 similarly to T* L in (16) with * rep¬ 
resenting u, y, or /, 

E id = diag (E id , E id , • • • , E id ), (39) 

L blocks 


Mr = 


M°r 




(40) 


and consists of the first L — r block-columns of 

M l 


Based on (36)-(38), we can write down the residual signal 
ffc l considering identification errors according to (17)- 
(19) and (27): 


f k,L =y k,L ~ T^y k,L - T^U ki L 

= T L A-r,L-r + ^,L + (r y L T") yfc.L 
+ (t£ - Tj?) u fc ,i 

= (f L)T - E id M x ) f, C _ T:i _ r + e fc>i ( 41 ) 


-M y r M“ 

y k,L 


u k,L 


true augmented fault signal f£_ T L _ T and the online I/O 
data z k,L- 

We regard f/_ T L _ T and z k,L as unknown but energy 
bounded. Hence f/_ T L and z k,L hr the first two terms 
of (43) lead to an estimation bias, while the online in¬ 
novation signal &k,L hr the third term causes zero mean, 
stochastic estimation errors. We would like to reduce 
the estimation bias by minimizing the matrix 2 -norms 
\\T. (G )\\ 2 i s — fi z )i an<4 at the same time minimize the 
Frobenius norm tr (GI e ,LG T ) by using the available in¬ 
novation covariance I e ,L- These three objectives are for¬ 
mulated by the following mixed-nornr problem: 

Gr off = arg min tr (t/E e rG T ) 

0 (44) 

s.t.E (Ta (G)T? (G)) < 1 2 S I, S = f,z 

where the matrix Q denotes the r-delay fault esti¬ 
mator (42), E denotes mathematical expectation over 
the identification innovations E; d , 7 / > 0 and i z > 0 
are the user-defined parameters to achieve a trade¬ 
off between estimation error variance and bias. Note 
that the matrix 2 -norms \\% (G)\\ 2 (s = f,z) are af¬ 
fected by the stochastic identification innovations E; d 
according to (43), hence their mathematical expecta¬ 
tions are used in (44). Note also that it is straightfor¬ 
ward to prove E (G) T s (G)) < 7 11 holds if and 
only if E (Ts(G)7?(G)) < ill in (44) holds. Here 
we use E (T s ( G ) 7? ( G )) in (44), because it brings a 
clear geometrical interpretation for parameter tuning 
as explained later in Section 5.2. With the tedious but 
straightforward derivations summarized in Appendix E, 
the above problem (44) can be explicitly written as 

Gr,oS = arg min tr (GI e .LG T ) (45a) 

0 


Similarly to Gn in (30), let the matrix Q denote the r- 
delay fault estimator based on the residual ffc.L, he., 

f(k - t) = Gr k , L - (42) 

It follows from (41) that the fault estimation error is 

A/(fc — T) = f(k — t) — In f ^k-T,L-T 

= — t7Ei d Mx — -Erif'j f k-T,L-T 

V_^_✓ 

T f (S) 

- £Ej d Mx z kfL + Qek,L 
TA0) 

(43) 

where I nf is defined in (29). It can be seen that E; d 
appears as multiplicative uncertainty coupled with the 


s.t. Q In, 


% -f L ,/ 


^ T 1 

-n ,r In, _ 


1 

H 


<7p 


(45b) 


GH z G t < ill, (45c) 

with Ilf and n, defined in (E.5) and (E. 6 ), respec¬ 
tively. The mixed-norm problem (45) can be easily trans¬ 
formed into an equivalent semi-definite programming 
(SDP) problem that can be solved efficiently [3]. Since 
the optimization problem (45) is determined only by the 
identification data and does not involve any online I/O 
data, it can be solved offline to obtain the robust fault 
estimator denoted as Gr,off- 


5.2 Parameter tuning using geometric interpretation 


Next, we will provide a systematic method to tune the 
two user-defined parameters 7 "j and 7 ^ by using a geo¬ 
metric interpretation of the mixed-norm problem (45). 
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With some matrix manipulations, we can see that the 
constraints (45b) and (45c) define two ellipsoids 


% 




(0-00)11/ (0-0o) T <0 O n / 0 o T -7 + 7 2j}, 

(46) 

ttz = {0|0LF0 T <7 2 z !} , (47) 


respectively, with 0o = X nj T J T Hj 1 . Since the objective 
function (45a) can be regarded as a measure of the dis¬ 
tance from Q to the origin 0 nf x(n -l), the optimization 
problem (45) is equivalent to finding the point 0 r , o ff in 
the set 17/ f) 17 z that is closest to the origin, as shown in 
Fig. 1. 


First, we would like to find the region of 7 “j and 7 ^ so that 
the optimization problem (45) is feasible and non-trivial. 
In the case that the origin 0 n/X ( n!j .L) € 17/OIL, we 
would have the trivial solution 0 r , o ff = 0 nf x(n y -L) which 
makes no sense for fault estimation. Hence 0 nf x(n y -L) ^ 
17 f and 11/^0 are both required, which implies the 
region of 7 / as below according to (46): 


1 \nin (0014/00 ) — 7/,min — 7/ < 1- (48) 


For a given 7 / satisfying (48), we solve the following 
optimization problem 

{0mini 7^ min } = &l'g min ^ 

<3,i l (49) 

s.t. (456) and (45c) 

whose solution gives the minimal 7 ^, referred to as 
7z min; that ensures 12/ P|17 z 7 ^ 0 . Therefore, we should 
select 7 1 G [7 Z m ini 00 ) to ensure feasibility of the op¬ 
timization problem (45). The ellipsoid 12- jln i n in Fig. 1 
represents the ellipsoid 12 - with 7 ^ = 7 ^ min , and its in¬ 
tersection with the ellipsoid 12 / includes only the single 
point 0 m in ■ 

By discarding the constraint (45c) from the problem (45) 
and fixing 7 / at the same given value as in (49), we 
formulate another problem 

0 i = arg min tr ( 0 S e>i 0 T ) 

Q (50) 

s.t. (456) 

Because the optimal solution 0i gives the shortest dis¬ 
tance from the origin to the ellipsoid 12 /, and more¬ 
over 0„ /X (n L) 44/, the solution 0i must lie at the 
boundary of the ellipsoid 12/, as shown in Fig. 1. Define 
7z,l — ^max (E (Tz (0i)7 z r (0i)))■ Let the ellipsoid 12 Z)i 
in Fig. 1 represent the set 12- with 7 ^ = 7 ^, and it has 
the solution Qi at its boundary. 


Similarly to the above obtained solution 0i of the 
problem (50), the solution 0 rjO ff of the problem (45) also 
lies at the boundary of the ellipsoid f2/. This allows the 
three terms of the fault estimation error in (43) to be 
explained using Fig. 1: 

1) The bias related to the first term 7/ (0) f k _ T L _ T is 
determined by the size of the ellipsoid 12 /; 

2) The bias related to the second term 77 (0) z kjL is 
determined by the size of the ellipsoid f l z (0 r ,off) 
with 0 r , o ff lying on its boundary, i.e., the ellipsoid 
n z with 7 f = Amax (E (T z (0r,off) T? (0r,off))); 

3) The fault estimation error variance related to the 
third term Qe k ,L is represented by the distance from 
the origin to the optimal solution 0 rjO ff- 

With the above basic geometric interpretation, we can 
analyze the performance trade-offs of the robust fault 

7/ jmin i l) and 72 G 
[7z mim 00 ) 1 as shown in Table 1. First, we fix jj and tune 
7 Z . In this case, the ellipsoid 17/ is fixed, which makes the 
first bias term in the first two rows of Table 1 remain con¬ 
stant. With the fixed 7 /, by increasing 7 ^ from 7 ? min to¬ 
wards 7 ^ 1 , the intersection set 12 / P| 12 - becomes larger, 
and the optimal solution 0 l jO ff moves from the point 0 m i n 
along the boundary of the ellipsoid 17/ towards the point 
01 . When we further increase 7 Z for 7 ); > 7 ^, the op¬ 
timal solution 0 r . o ff of the problem (45) would remain 
located at the point 0 i, because Q\ satisfies both con¬ 
straints (45b) and (45c) and gives the shortest distance 
to the origin according to the problem (50). Therefore, 
the size of the ellipsoid 17- (0 r ,off), which determines the 
second estimation bias term in the first two rows of Table 
1 , monotonically increases for 72 G [ 7 ^ min , 7 ? : ) and re¬ 
mains constant for 72 G [ 7 ^, 00 ). The distance from 
the origin to Gr.offi which determines the fault estima¬ 
tion error variance in the first two rows of Table 1, mono¬ 
tonically decreases for 7 ^ G [7? ;ln i n , 7 Z) i) and remains 
constant for 7 ^ G [ 7 Z 1 , 00 ) . For the third row of Table 
1 , we tune 7 / and select a sufficiently large value of 7 Z 
that ensures the problem (45) to be feasible. With 7 / in¬ 
creasing, the size of the ellipsoid 17/, which determines 
the first bias term in the third row of Table 1, monoton¬ 
ically increases. Meanwhile, the optimal solution 0 r , o ffi 
which lies at the boundary of the ellipsoid 17/, moves 
closer to the origin. Therefore, both the second bias term 
and the fault estimation error variance in the third row 
of Table 1, which are determined by the size of the ellip¬ 
soid 12 z (0r jO ff) and the distance from the origin to the 
point 0r, o ff, monotonically decrease. 

We summarize the data-driven robust design in Algo¬ 
rithm 2. The nominal design Q n obtained from Algo¬ 
rithm 1 can be used as a benchmark for tuning 71 and 72 
in Step 2 of Algorithm 2. For example, compared to the 
nominal design, the robust design achieves smaller aver- 


estimator 0 riO ff when tuning 7 / G 



Table 1 

Trade-offs between fault estimation bias and error variance of the robust fault estimator C/ r ,off at time instant k when tuning user- 
defined parameters 7 ^ and 7 ? in (45): “Constant”, and means that the performance criterion in the corresponding 
column remains constant, monotonically increases, and monotonically decreases with regard to the user-defined parameter 
specified in the corresponding row, respectively. 


User-defined 

parameters 

First bias term 

ETHMfL/-! 

Second bias term 

E \\Tz (Gr,o ff) Zk.L || 2 

Variance 

tr (Gr,oS^e,LGj' 0 g) 

ll € [7^,min, t2,i] 

Constant 

/ 

\ 

ll € [ 7 ? 1 , 00 ) 

Constant 

Constant 

Constant 

7/ G [7/,min; l) 


\ 

\ 



6 Data-driven robust receding horizon fault es¬ 
timation with online optimization 

The online I/O data is regarded as unknown in Algo¬ 
rithm 2. In order to better exploit the available online 
data, this section proposes an online mixed-norm opti¬ 
mization approach. This can further reduce the estima¬ 
tion errors when the online I/O data have large ampli¬ 
tudes, at the expense of increased computational burden. 

6.1 Online mixed-norm problem 


Fig. 1. Geometric interpretation of the mixed-norm problem 
(45): the constraints (45b) and (45c) define the ellipsoid f2/ 
centered at Go and the ellipsoid fl z centered at the origin O, 
respectively. Lying at the boundary of the ellipsoid SI/, the 
optimal solution Gr,o ff gives the shortest distance measured 
by the objective function (45a) from the origin to the inter¬ 
section set SI/ With 7 ? = 7 ? ;m i n , the ellipsoid fl z be¬ 

comes SIz im i n (S/ m in) in green which intersects with the ellip¬ 
soid SI/ at a single point Gmin- At the boundary of the ellip¬ 
soid SI/, Q\ gives the shortest distance from the origin to the 
ellipsoid SI/. The ellipsoids ff Zi i (Gi) in blue and SI Z (<5 r ,off) 
in red represent the ellipsoids SI Z with Gi and Gi,oS lying at 
the boundary, respectively. 


With the notation 

Pk,L = M/z/.^, 

we divide 3 k / into L row blocks as in 


Pk,L = 


n T 


Ph ••• Pk,L 


(51) 


(52) 


with £ R w . Then the term C/EjdM^z^/, in (43) can 
be rewritten as 


aged worst-case bias if 7 ^ < A max (E ( T s ( G n ) T7 (£n))) 
(s = /, 20- 


Algorithm 2 Data-driven robust RH fault estimation 

1) Complete the steps 1-3 in Algorithm 1 ; compute 

Mf. and M[ according to (13) and (35). 

2) Tune £ 7/,mho 1 ) and 7 1 G [7z,min>°°) ac¬ 
cording to the performance trade-offs shown in 
Table 1, where 7 1 min and 7 / min are obtained from 
the optimization problems (48) and (49) respec¬ 
tively. 

3) Solve the problem (45) to compute the robust RH 
fault estimator Q T G ff. 


t/E id M2Z fcji = G^idPk,L 


Ei d/ dfcp 


Pk,l ^ Ifly 

Ei d /3fe,2 

= G 

Pk ,2 ® Iriy 

. E id /3fc,i _ 


.Pk,L ® Iriy . 


Fk,L 


(53) 


according to the property of Kronecker product [4], 
Based on (53), the estimation error in (43) becomes 

A f(k — t) =Tf (G) f/_ T l _ t - GT k ,Lve c (E id ) + Ge k , L - 

(54) 

Then the statistics of vec (Ej d ), i.e., 

E (vec (E id ) vec (E id ) T ) = I N ® E e , 
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can be exploited to evaluate the fault estimation error 
variance. Therefore, we formulate the following opti¬ 
mization problem similarly to (44): 

<?r,on = arg min tr (gE e ^ L G T + GT k}L (I N 0 E e ) 
s.t. e(T/ (0)77(0)) < 7 ^ 

. (55) 

with the user-defined parameter jf. The constraint in 
the above optimization problem (55) can be explicitly 
written as (45b). The optimization problem (55) has to 
be solved at each time instant to update the robust fault 
estimator t/ r;0n because T kt L in the cost function is de¬ 
termined by the online I/O data according to (51)-(53). 

6.2 Parameter tuning using geometric interpretation 

Since the online mixed-norm problem (55) has the struc¬ 
ture similar to that of the offline mixed-norm problem 
(45), the performance trade-offs by tuning 7 f in (55) are 
also similar to that explained in Table 1. 

The proposed data-driven robust fault estimation with 
online optimization is summarized in Algorithm 3. In 
order to reduce the computational burden of online op¬ 
timization, the problem (55) is implemented only if the 
estimation bias of the offline designed fault estimator is 
larger than a user-defined threshold a, as shown in Step 
2 of Algorithm 3. 


studied in [12-14,16]: 

4-c(l) — ^cXcit') T -H c U c (t), 

y c {t) = C c {t ), 

r-0.0366 0.0271 0.0188 -0.45551 
A — 0.0482 -1.01 0.0024 -4.0208 

0.1002 0.3681 -0.707 1.42 

L 0 0 1 0 J 



r 0.4422 0.1761 1 

i 

'1 0 0 0' 

to 

0 

II 

3.5446 -7.5922 
-5.52 4.49 

ii 

0 10 0 

0 0 10 


. 0 0 

i 

0 111 


With a sampling rate of 0.5 seconds, the discrete-time 
model (2) is obtained, with D = 0 and F = 1 4 . The 
process and measurement noises, w(k) and v(k), are as¬ 
sumed to be zero mean white noises, respectively with 
covariances of Q = 0.1 6/4 and R = O. 64 / 4 . 

Since the open-loop plant is unstable, an empirical sta¬ 
bilizing output feedback controller is used [ 12 ], i.e., 

u( k ) = -[00 -on - 0.1 ] y( k ) + v( k ), ( 56 ) 

where r](k) is the reference signal. 

In the identification experiment, the reference signal 
r)(k) is zero-mean white noise with the covariance of 
diag (1,1), which ensures persistent excitation. We col¬ 
lect N = 1000 data samples from the identification 
experiment. In the identification algorithm, the past 
horizon is selected as p = 10 . 


The offline designed fault estimator Gr,ott from Algo¬ 
rithm 2 can be used as a benchmark for tuning in 
Step 2.2 of Algorithm 3. For example, compared to G T ,oS, 
the online optimization (55) achieves smaller averaged 

worst-case bias if < A max (e (Tf (<? r , 0 ff) Tj (<? r ,off)))- 


Algorithm 3 Data-driven robust RH fault estimation 
with online optimization 


1) 

2 ) 


Follow Algorithm 2 to compute the offline de¬ 
signed fault estimator G T>0 g. 


If A min (E (T z (Sr.off) 77 (Sr,off))) ||z fc>L ||* > « (<* 
is a user-defined threshold), the online optimiza¬ 
tion in the following steps is implemented; other¬ 
wise, the offline designed estimator G r ^ 0 g is used. 
2.1) Compute T k ,L according to (51)-(53). 


2.2) Tune 7 “j 


7/,min* 1 


similarly to Step 2 of 


Algorithm 2, with 7 ^ in defined in (48). 

2.3) Solve the problem (55) to compute the robust 
RH fault estimator t/ r , 0 n- 


The considered fault cases include: 

• Actuator faults: E = B, G = D, 

• Sensor faults: E = 04 X 2, G = [J 1 0 o] T - 

The case of simultaneous actuator and sensor faults is 
not included here, because all the considered algorithms 
can be applied to the simultaneous scenario in a straight¬ 
forward way, and their performance comparisons are the 
same as in the case of separate actuator or sensor faults. 

The simulated fault signals in both fault cases are the 
same: 


m 



sin ( 0 .l 7 rfe) 



0 < k < 50, 
k > 50. 


We will compare the following fault estimation methods: 


7 Simulation studies 

Consider a continuous-time linearized vertical take- 
off and landing (VTOL) aircraft model that has been 


• AlgO: the RH fault estimator using accurate Markov 
parameters, described in Section 4. 

• DONG: the method proposed by [12]. 

• Algl: the data-driven nominal RH fault estimator Q n 
proposed in Algorithm 1 ; 
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Actuator fault 


Sensor fault 


• Alg2: the data-driven robust RH fault estimator Q It0 s 
proposed in Algorithm 2; in Step 3 of Algorithm 2, we 

select 7 ^ = A max (e (Tj (Qn) Tj (Sn)'j'j , and 

7z = 0-5 (7z, min + 7*,i) • (57) 

• Alg3: the data-driven robust RH fault estimator 5 r ,on 
with online optimization, proposed in Algorithm 3; 
in Step 2 of Algorithm 3, we select a = 300 as the 
threshold to determine whether or not the online opti¬ 
mization should be implemented; is set to the same 
value as in Alg2. 

We select the estimation horizon length L = 30 for the 
considered five algorithms. 


True fault 

- AlgO 

DONG 

- Algl 

-Alg2 

- Alg3 


-4; 





150 200 250 300 

time instant (k) 


In order to show the necessity of compensating for the 
identification errors, we make the identification-error- 
effect term 77 (Q) ■ z^l in (43) significantly large by set¬ 
ting T](k) = 15. Fault estimates from the above five al¬ 
gorithms are illustrated in Fig. 2, and the distributions 
of their fault estimation errors are shown in Fig. 3. By 
using accurate Markov parameters, AlgO achieves unbi¬ 
ased fault estimation in both fault scenarios. Note that 
DONG cannot be directly applied to sensor faults in the 
unstable open-loop VTOL model [12], hence it is not in¬ 
cluded in Fig. 2 and 3(b) for sensor faults. Because of 
neglecting the effect of identification errors, both Algl 
and DONG yield estimation biases even larger than the 
amplitude of true faults. In comparison, Alg2 obtains 
its robustness to identification error by solving an off¬ 
line mixed-norm problem, as shown in Fig. 3(a). How¬ 
ever, the poor performance of Alg2 in our sensor fault 
case (Fig. 3(b)) shows the limitation of neglecting the 
online availability of I/O data in the offline mixed-norm 
problem. Compared to Alg2, Alg3 significantly reduces 
estimation bias, as shown in Fig. 3(b), by formulating 
an online mixed-norm problem to exploit online I/O 
data. This performance improvement is achieved at the 
cost of higher online computational burden. When im¬ 
plemented with YALMIP [21] in the MATLAB2011b 
environment, on a computer with a 3.4 GHz processor 
and 8 GB RAM, the averaged and peak computational 
time per sample of Alg3 are 1.70s and 2.05s for the es¬ 
timation horizon length L = 30, while those of Alg2 are 
8.37 x 10 _6 s and 3.17 x 10 _5 s respectively. We will inves¬ 
tigate the computational efficiency of Alg3 for real-time 
implementation in future work. 

In order to illustrate the performance trade-offs of 
Alg2, we set 7 ? as in (57) and tune 7 ^ under the con¬ 
dition of different reference signals ij(k). Fig. 4 shows 
how the fault estimation bias, error variance and root 
mean square error (RMSE) vary with 7 ‘j, which can 
be explained as follows using Table 1. According to 
the fault estimation error analysis in (43), the fault 

estimation bias is related to both 7/ (f? r ,off) f(L r l-t 


Fig. 2. True fault signal and fault estimates from different 
algorithms. 


AlgO 

DONG 



.3 1 -.-.-.-■-■-■-■-' 

-2 -1 0 1 2 3 4 5 6 

AW - /iW 

(a) Actuator faults 



-35-■-■-■-■ 

-50 0 50 100 150 

AW-/1W 

(b) Sensor faults 

Fig. 3. Distribution of fault estimation errors when 
r/(k) = 15. Circles: 1000 estimation errors by dif¬ 

ferent fault estimation algorithms based on 1000 on¬ 
line I/O data samples. Ellipses: the 3<r-contour of 
the approximated two-dimensional Gaussian distribution 
of the 1000 estimation errors, i.e., the contour at 

[/(*0 - /(*=)] cov" 1 (/(fc)) [/(fc) - /(&)] = 3. 

and 77 (Sr,off) For r](k) = 0 or rj(k) = 1, the 

online I/O data z^.l have small amplitude, thus the 
total estimation bias is dominated by the bias related 
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to 7/ (Gr.off) l _ t which monotonically increases 
with 7 ‘j according to the third row of Table 1. This 
explains the fault estimation bias curves for ry(fc) = 0 
and r](k) = 1 in Fig. 4. For tj(k) = 2, the online I/O 
data z k t L have relatively large amplitudes, hence for 
relatively small values of 7 ^ the total estimation bias is 
dominated by the bias related to T z (l?r,off) z k,L which 
monotonically decreases with 7 ‘j, and for relatively large 
values of 7 ^ the total estimation bias is dominated by 
the bias related to 7/ (G r ,off)f i- r ,L- T which monoton¬ 
ically increases with 7 ^, according to the third row of 
Table 1. This explains the fault estimation bias curve 
for r](k) = 2 in Fig. 4. The nronotonic decrease of the 
fault estimation error variances with yj can be directly 
explained with the third row of Table 1. As the objec¬ 
tive function of the optimization problem (45), the fault 

estimation error variance tr ^<?r,offS e ,L(?r I off) for dif¬ 
ferent reference signals ry(fc) is the same because it does 
not depend on the reference signal r/(k). Combining the 
increase of fault estimation bias and the decrease of 
fault estimation error variance with y“j, there exist the 

optimal 7 y * € ^ 7 y min , 1^ such that the RMSE achieves 
its minimal value, as can be seen in Fig. 4. It is also 
shown that the minimal RMSE is achieved at a larger 
value of y“j * when the amplitude of rj(k) increases, 
because the online I/O data have larger amplitudes 
with larger rj(k), thus the decrease of the bias related 
to T z (t/r,off) z k.L dominates the fault estimation bias. 
Based on the above insights, we can anticipate how the 
estimation performance of Alg2 varies with different y% 
for a fixed 7 ^, as well as the performance trade-offs of 
Alg3. Their performance curves are not plotted due to 
the space limitation. 



Fig. 4. Estimation performance of Alg2 when tuning ■y'j under 
different reference signal r/(k) 

be applied to sensor faults of an unstable open-loop plant 
which could not be directly addressed previously. In the 
formulated RH fault estimator, the identification errors 
appear as multiplicative model uncertainty coupled with 
the unknown faults and the online I/O data. Then, two 
mixed-norm problems were formulated to enhance ro¬ 
bustness. One can be solved offline by regarding the on¬ 
line I/O data as unknown signals. The other further re¬ 
duces estimation errors for larger I/O data by exploiting 
their online availability in the mixed-norm problem, and 
it requires online optimization. Based on geometric in¬ 
terpretations of the mixed-norm problems, systematic 
methods were given to tune the user-defined parameters 
therein. Comparisons using a simulated aircraft model 
illustrated the advantages and the effectiveness of our 
proposed method. 
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8 Conclusions 

This paper has investigated data-driven fault estimation 
and its robustness against stochastic identification er¬ 
rors. First, we proposed an RH fault estimator that can 
be parameterized with the predictor Markov parame¬ 
ters. Its condition for unbiasedness generalizes that of a 
recently reported data-driven fault estimation method. 
An immediate benefit is that our proposed method can 
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A Lemmas for Theorem 1 


Lemma 1 Define x e (0) € R”, f e (i ) € R n/ , and r e (i) G 
R"» (i > 0) as the initial state, input and output signal 
of the fault subsystem (&,E,C,G), respectively. There 
exists a non-zero initial state £ e (0) such that r e (0 ) = 
r e (l) = • • • = r e (L) = 0 for all L > v + t, if and only if 


(i) O r x e ( 0) = 0; 

(ii) the system 


fik + l) = \$ - E (H() C$ T ]x e (k) 


K d 


(A.l) 


r e (k) = I -H( (Hi) C$ T x e (k) 


is unobservable; 


13 



(in) the inputs {f e (i)} take the form 


B Proof of Theorem 1 


fe(i) = ~ (H() C^ T K l d x e { 0). (A.2) A solution f£_ T L _ T to the problem (20) satisfies 


In Lemma 1, r e (0) = • • • = r e (r — 1) = 0 is ensured 
because of the condition (i) and the zero Markov ma¬ 
trices according to Assumption 2, 

while r e (t) = • • • = r e (L) = 0 is ensured by the condi¬ 
tions (ii) and (iii). Lemma 1 can be proved by slightly 
modifying Lemmas A.l and A.2 in [20]. 


*I t Z- 1 l *L,t%- t ,L-t = *L ,rS"ir fcii . (B.l) 

Let Af£_ T)i _ T = ff_ Tii _ r - ffc_ TjL _ T denote the esti¬ 
mation error. By substituting (19) into (B.l), we have 


From Lemma 1 we can see that perfect reconstruction 
of system inputs {f e (i)} from system outputs {r e (z)} is 
impossible if the unobservable input signal (A.2) is non¬ 
zero. Hence, next, we will investigate the link between 
the unobservable input signal (A.2) and the system prop¬ 
erty of ($, E, C, G). 

By setting i = 0, (A.2) becomes 


fe (0) = - {Hi) C<S> T x e ( 0). (A.3) 

Then, according to the condition (i) and the unobserv¬ 
ability of the system (A.l), there must exist a scalar A 
and a non-zero x e (0) such that [37] 


K,i-XI 

o T 


X e (0) 


•S>-XI E 
Or 0 


*e(0) 
fe (0) 


' <S>-XI E 
. Or + l H( 

where H l defined in (24) equals to 


» e (o) 

f e (0) 


= 0 , 

(A.4) 


o 

Hi 


because 


H t { , H{ , • • ■ ,Hl_ 1 are zero matrices according to As¬ 
sumption 2. With (A.3) and (Kd — A/)x e (0) = 0 in 
(A.4), we can rewrite f e {i ) in (A.2) as 


* 1 , 




Af; 


k—T,L—' 


= 'fTs -1 


L e k,L, 


which implies ^I )T S e ^ L>T E (Af£_ TL _ T ) = 0 by 

taking expectations on both sides. Therefore, the un¬ 
biasedness condition of the estimate in (23) reduces to 
the analysis of the linear equation 

*l,t E (Af *_ rtL _ r ) = o (B.2) 

since TV =A f{^ L ,r)- 

The rest of the proof follows the intuitive arguments 
below. According to Lemma 1, (A.5), and the definition 
of ff_ T l _ t in (19), there are three scenarios: 


1) When ($, E, O t+ i, H;f) has no invariant zeros, 
the non-zero initial state x e (0) in Lemma 1 does 
not exist according to (A.4), thus (B.2) implies 

E ^Af)?_ TjL _ r ^ = 0, i.e., unbiased fault estimation. 

2) When ($, E, O t+ \, H^) has invariant zeros, (B.2) 
implies that for each invariant zero A, the expected 
error of the r-delay fault estimate f{k — r) is 


f e (i) = A7e(0). (A.5) 


E (A/(fc — r)) = A L_r_1 E (A/(fco)) (B.3) 


The above analysis indicates that the unobservable in¬ 
puts {f e (i) = A*/ e (0)} are determined by the invariant 
zero A of (<&, E, O t +i, H(), as shown in the following 
lemma: 

Lemma 2 Considering the non-zero initial state x e (0) 
in Lemma 1, there are two types of the invariant zeros 
A of the fault subsystem (i>, E, O t +i, H£) in (A.4): 1) A 
is an unobservable mode , then (A.4) implies f e ( 0) = 0, 
thus the input signal {f e (i) = A*/e(0)} is constantly zero; 
2) X is a transmission zero, then f e { 0) ^ 0, thus the 
unobservable input signal {f e (i) = A*/e(0)} is non-zero. 

Lemma 2 directly extends Lemmas 1 and 2 in [35] which 
considers only the case r = 0. 


in the estimation horizon [ko, k] (ko = k — L + 1). 

2.1) If all the invariant zeros of (<f>, E,O r + i,H£) 
correspond to unobservable modes, it follows 
from the case 1) in Lemma 2 that the ex¬ 
pected estimation error (B.3) is zero because 
E(A/(fc o )) = 0. 

2.2) If transmission zeros exist but are all stable, 
i.e., |A| < 1, it follows from the case 2) in 
Lemma 2 that E(A/(fco)) 7 ^ 0 and the ex¬ 
pected estimation error (B.3) asymptotically 
reduced to zero as L goes to infinity. 

The scenarios 1) and 2.1) correspond to the case (i) of 
Theorem 1, and the scenario 2.2) corresponds to the case 
(ii) of Theorem 1. 
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C Proof of Theorem 2 


equivalent: 


For the original system model (2), the extended output 
equation in the time window [fco, k] is 

y k,L = 0LX(k O ) + k.L + fk,L + + v fc,Zo 

(C.l) 

where 0 L , &[, and ^ are defined in the same 

way as Ol and T£ in (16). According to (C.l), we can 
rewrite (17) and (18) as 

r k,L = {I ~ T|) (y fcjL - 3 ?^ Ufc.L) 

= (/ - T" ) (Vi^fco) + &lh,L + ^Z°w fciL + v* >L ) 
= (J-T»)[^ - ' er.,.. 



by following the relation between the original system 
model (2) and its predictor form (3). Similarly to T J L 
in (19), ST[ in (C.2) consists of the Hrst L — r block- 
columns of . 

Define f fc;L = y fe;L - u k%L and 


t L = cov (^ w w k>L + v k>L ). 

Comparing (19) with (C.2) leads to 

r k ,L = (I-T y L )r k , L , * L ,r = (I- TD'f'L.r, 

Ze,L = (I-T y L )£ L (I-Tl) T . 

Then by substituting (C.3) into (22), the estimate of 
f k-r,L-r becomes 


= (*L 


S7 


which is actually the LS estimate proposed in [35] based 
on the original system model (2). 


D Proof of Theorem 3 

Split T j l t into two blocks as T^ r t{ t , with T^ t 
consisting of the first L — t — 1 block-columns of T{ , 
and T£ consisting of the last block-column of T{ . 
With these notations, unbiased fault estimation can be 
proved by showing that t£ t E(A f(k — r)) = 0 because 
T^ has full column rank according to Assumption 2. 

According to (26), the following two expressions are 


e€K([o L t{ iT ])n^(T{, r ), (D.l) 

^([h^ (D.2) 

Since the two sufficient conditions for (asymptotic) unbi¬ 
asedness in Theorem 1 imply e = 0 and e —> 0 (L —> oo) 
for (D.l), it then follows from the equivalence between 
(D.l) and (D.2) that the sufficient conditions in The¬ 
orem 1 also imply e = 0 and £ —> 0 (L oo) for (D.2), 

or equivalently, 1Z (t{ = {0} and 1Z ^t£ r ^j — t {0} 

(L —> oo). Therefore we can conclude that the sufficient 
conditions in Theorem 1 imply (asymptotically) unbi¬ 
ased fault estimation for (D.2). Similarly, we can prove 
the necessary condition for the (asymptotically) unbi¬ 
ased fault estimation. 


E Computation of E (T s (G) Tj (G)) 

By dividing Mx in (40) into L row blocks as 

Mr = T , (E.l) 

With M T) j G ]Rn/x(m.n u +(L-T)n / ) ) we de fi ne p x as 

tr^Mr.iM^i) ••• tr^MT.iM^ 

= 

tr(M T i) tr(M T lM? 2 ) — tr(M T z,M^ l ) 

’ (E.2) 

P, is defined similarly to (E.2), by dividing in (41) 
into L row blocks as in (E.l). Then, 


E (T, (S) 7} T (S)) = [ g x„, ] [ J L ’ T ] f 

L L,r -Lrif J LAif 

(E.3) 

e(t z (G)t? (G)) =gn z g T (e.4 ) 

with 

n/ = f L ,rn, T + E (E id M T M?ET) 

= f iiT fI iT + P x ®E e , 

n, =E(E id M7(M|) T ET) =P z g,E e . (E.6) 
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